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Abstract: Feedback circuits in biochemical networks which underly cellular sig- 
naling pathways are important elements in creating complex behavior. A specific 
aspect thereof is how stability of equilibrium points depends on model parameters. 
For biochemical networks, which are modelled using many parameters, it is typ- 
ically very difficult to estimate the influence of parameters on stability. Finding 
parameters which result in a change in stability is a key step for a meaningful 
bifurcation analysis. We describe a method based on well known approaches from 
control theory, which can locate parameters leading to a change in stability. The 
method considers a feedback circuit in the biochemical network and relates stability 
properties to the control system obtained by loop-breaking. The method is applied 
to a model of a MAPK cascade as an illustrative example. 
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1. INTRODUCTION 

Feedback circuits are an important structural fea- 
ture of biochemical networks [Tyson and Othmer, 
1978]. The presence of complex behavior such as 
bistability, i.e. the existence of several stable equi- 
libria, and sustained oscillations can be attributed 
to the presence of feedback circuits [Cinquin and 
Demongeot, 2002, Kaufman and Thomas, 2003]. 
These types of complex behavior are directly re- 
lated to how feedback circuits influence stability 
properties of equilibria. In consequence, stability 
analysis of biochemical networks involving feed- 
back is a recurring field of interest, and several 
theoretical results have been obtained [Dibrov 
et al, 1982, Thron, 1991, Angeli and Sontag, 
2004]. 

Models for biochemical networks in cellular signal- 
ing typically contain a large number of parameters 



whose values are not exactly known and which 
can even vary due to differential gene expression 
(e.g. the concentration of an enzyme) or external 
influences (e.g. cofactors). These parameters often 
have a considerable influence on stability, which 
needs to be evaluated in order to understand the 
function of a network [Eifiing et al., 2007, Kim 
et al., 2006]. 

A classical tool to study the influence of parame- 
ter variations on stability is bifurcation analysis. 
It has been applied to many cellular signalling 
systems, such as the lac operon [Yildirim and 
Mackey, 2003] and the MAPK cascade [Marke- 
vich et al., 2004, Chickarmane et al., 2007], to 
name but a few. When considering models with 
many parameters, one faces the difficulty that in 
classical bifurcation analysis, only one parameter 
at a time can be varied. Thus the effect of simul- 
taneous variations in several parameters can not 
be evaluated properly. 
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In this paper, we present a new approach to lo- 
cate bifurcations in systems with feedback loops 
containing many parameters which may be varied 
simultaneously. To this end, we make use of an ap- 
propriate frequency domain description of the sys- 
tem and of mathematical conditions representing 
necessary conditions for a bifurcation. The paper 
is structured as follows. In Section 2, we present 
the theoretical results required for our approach 
and suggest an optimization-based method to ac- 
tually find interesting parameter values. In Sec- 
tion 3, we apply these results to a model of the 
MAPK cascade with a negative feedback circuit 
[Kholodenko, 2000]. The relevance of our results 
is discussed in Section 4. The mathematical proofs 
of the theoretical results are not presented in this 
paper but will be provided elsewhere (Waldhcrr 
and Allgower, in preparation). 



2. THEORETICAL BACKGROUND 

2.1 The loop-breaking approach 

We consider a nonlinear differential equation 
which may describe the biochemical network con- 
stituting a cellular signaling pathway 



£ : x = F(x,p), 



(1) 



with x G M™ and p G V, where V is a con- 
nected subset of R m . Typically x will represent 
the concentrations of the signaling molecules, and 
p collects parameters like reaction constants or 
enzyme concentrations. We assume that an equi- 
librium point x(p) exists for all parameter values 
and can be computed at least numerically, such 
that F{x{p),p) = for all p G M m . 

Mathematically the system (1) is said to contain 
a feedback circuit if the influence graph of its 
Jacobian ^ contains a nontrivial loop [Cinquin 
and Dcmongcot, 2002]. We want to study the in- 
fluence of such a feedback circuit on the dynamical 
properties of the system. Control theory provides 
efficient tools to study this problem. A useful 
approach in our setup is to consider the system 
(1) as a closed loop control system. It is then 
possible to study the corresponding open loop 
system, and one can resort to the rich stability 
theory developped for control systems. 

An open loop control system corresponding to 
the closed loop system (1) is obtained by loop 
breaking, as defined in the following. 

Definition 1. A loop breaking for the system (1) 
is a pair (/, h), where / : 1" x I x R m -> R™ is a 
smooth vector field and h : R n — > R is a smooth 
function, such that 



(3) 



The corresponding open loop system is then given 
by the equation 

x = f(x,u,p) 
y = h(x). 

The closed loop system can again be obtained 
by "closing the loop", i.e. setting u = y. Notice 
that by the assumption that an equilibrium exists 
for the closed loop system, the open loop system 
also has the equilibrium x(p) when choosing u = 
h(x(p)). This input is denoted as u(p) — h(x(p)). 

Since our main interest is in stability properties 
of the equilibrium point x(p), we can restrict the 
analysis to the linear approximation of the sys- 
tems (1) and (3) around the equilibrium point. By 
using Laplace transformation, the linear approxi- 
mation of the open loop system (3) can be repre- 
sented by a linear parameter-dependent transfer 
function 



G{p,s) = 



(4) 



r(p, s) 

where q and r are polynomials in the complex 
variable s with coefficients depending on p. As 
a technical restriction, we assume that the open 
loop system has no poles or zeros on the imaginary 
axis, i.e. r(p, •) and q(p, •) are assumed to have no 
roots on the imaginary axis for any value of p G V 
throughout this section. 

2.2 Properties of the closed and open loop systems 

Stability of an equilibrium point of the closed loop 
system depends on the position of the eigenvalues 
of the Jacobian ^(x(p),p). To characterize these 
eigenvalues from conditions on the open loop 
system, we have the following theorem. 

Theorem 1. Let A(p) = || (x(p),u(p),p) and 
A c i(p) = ^{x{p)tP). Assume that s G C is not 
an eigenvalue of A(p). Then s is an eigenvalue of 
A c i(p), if and only if G(p, s Q ) = 1. 

The proof of Theorem 1 is based on a representa- 
tion of G as 

det(sl - A cl (p)) 



G(p, s) 



det(sJ - A{p)) 



+ 1. 



F{x,p) = f{x,h{x),p). 



(2) 



Parameter values on the border of stability are 
characterised by the matrix A c i{p) having eigen- 
values on the imaginary axis. To study the corre- 
sponding property in the frequency domain rep- 
resentation of the open loop system, we introduce 
the notation of critical frequencies and gains. 

Definition 2. We say that io c G R is a critical 
frequency and k c G R a corresponding critical gain 
for the transfer function G(p, •) (4), if 

k c q(p,ju, 



1. 



(5) 



In general, different critical frequencies and gains 
will be obtained for different values of p. 

The critical frequencies can be characterized inde- 
pendently of the critical gains. This result follows 
from (5), because the transfer function value at 
the critical frequency has to be a real number. 

Proposition 1. lo c is a critical frequency for G(p, •), 
if and only if 

lm(q(p, juj c )r(p, -ju c )) = 0. (6) 

There exists a unique corresponding critical gain 
for any critical frequency lo c , which is given by 



for all parameters p G V. Thus, it can be used to 
characterise a set of critical frequencies as being 
minimal. 

Definition 3. Under the assumptions of Prop. 2, 
the set of critical frequencies fi c (p) is called mini- 
mal, if it contains exactly the minimal number of 
elements according to Prop. 2. 

If fl c (p) is minimal, we can label the roots of 
(6) in a consistent way, and write O c (p) = 
{u>l(p),u)c(p), . . . , u)"(p)\, where the lo 1 c can be 
identified with different solution branches of the 
polynomial equation (6). 



The equation (6) is a polynomial in w c , thus all 
critical frequencies can be computed numerically 
for fixed parameters p. The set of all critical 
frequencies for the transfer function G(p, •) is 
given by 

n c (p) = {to G R | lm(q(p,jui)r(p, -jw)) = 0} 

(8) 

The concept of critical frequencies and critical 
gains can be understood intuitively when consid- 
ering the Nyquist plot of the transfer function 
G{p 1 juj). A critical frequency is any value u> at 
which the Nyquist plot crosses the real axis. The 
corresponding critical gain is the value k(p) = k c 
which scales the Nyquist plot in such a way that 
the crossing point at the critical frequency is 
mapped to 1 in the complex plane. However, this 
intuitive way of scaling the Nyquist plot would 
require to keep all critical frequencies in f2 c (p) 
constant when varying parameters, which would 
be a strong restriction. The next section presents 
an approach to overcome this restriction. 

2.3 A minimal set of critical frequencies 

The number of critical frequencies that exist for 
a given open loop system is often predefined by 
the position of the open loop poles and zeros in 
the left or right half complex plane. The following 
proposition guarantees the existence of a minimal 
number of critical frequencies. 

Proposition 2. Let a = \p + —p-+z-—z + \, where 
p + (p-) is the number of poles of G(p, •) in the 
right (left) half complex plane and z + (z_) is the 
number of zeros of G(p, •) in the right (left) half 
complex plane. Then, for any p G V, £l c (p) has at 
least a distinct elements, if a is odd, and at least 
a — 1 distinct elements, if a is even. 

Since we assumed that G{p, •) has no poles or zeros 
on the imaginary axis, the number a is the same 



2.4 Existence of critical parameter values 

The concept of critical frequencies and gains is 
now applied to the problem of how stability de- 
pends on parameters. We study the problem of 
finding critical parameters p c G V on the border 
of stability, i.e. such that the eigenvalues of the Ja- 
cobian A c i(p c ) are located on the imaginary axis. 
Then there exist typically parameters po and p\ 
in a neighborhood of p c such that the equilibrium 
x{pi) is stable and x{p2) is unstable. 

The following theorem uses the loop-breaking 
approach and the concept of critical frequencies to 
characterise the existence of critical parameters. 

Theorem 2. Assume that f2 c (p) is minimal for all 
p G V . Then there exists p c G V such that jui l c (p c ) 
is an eigenvalue of A c i(p c ), if and only if there 
exist po,pi G V such that G(po, ju) l c (po)) < 1 and 
G(pi,ju l c (pi)) > 1, for w'(-) G Q c (-) and for any 
i G {1,2,..., a}. 

Thus, instead of having to look at how the n 
eigenvalues of the closed-loop system change with 
parameters, we have reduced this to one number, 
given by G{p, ju % c (p)), which contains all informa- 
tion about whether the system changes its sta- 
bility properties when changing parameters. The 
result is global in the sense that the parameters 
Po and pi can be arbitrarily far apart from each 
other, still under the given conditions existence of 
critical parameters p c is guaranteed. 



2. 5 Searching for critical parameter values 

The theoretical approach outlined above can be 
used to search for parameter values such that 
the equilibrium point x(p) of the closed loop 
system (1) changes its stability. For a biochemical 
system, there are often nominal parameters po, 
giving rise to the equilibrium point x(po). We want 



to find parameters p\ such that x(po) and x{p\) 
have different stability properties. 

In view of the methods presented in this paper, 
given the open loop transfer function (4) , one first 
needs to identify the critical frequency that is to 
be considered. This choice depends on the type of 
stability change one is looking for. When taking 
uo c = 0, it is possible to search for zero eigenvalues, 
and if u> c > 0, nonzero imaginary eigenvalues may 
be encountered, typically giving rise to a Hopf 
bifurcation in the closed loop system. The nominal 
transfer function value at the critical frequency 
is G(po,ju> c (po))- To change stability properties, 
we will then define a value 7 as either 7 > 1, if 
G(po, juj c (po)) < 1, or as 7 < 1 otherwise. Then, 
any solution to the nonlinear equation 

G(pi,iw c (pi))=7 (9) 
gives parameters p\ such that x(po) and x{p\) 
have different stability properties as indicated by 
the chosen critical frequency lo c . This method has 
been implemented using a nonlinear constrained 
optimization algorithm from the Matlab Opti- 
mization Toolbox [The MathWorks Inc., 2006]. It 
allows to efficiently compute parameter values for 
the desired transfer function value for medium 
sized systems, as the example presented in the 
following section illustrates. 

Once a parameter p\ is known, we can use a 
straight line going from p to pi, defined as 
Pa = Pa + Mpi — Po)- The change in dynamical 
behaviour along this line can then be studied 
using classical bifurcation analysis with respect to 
/j, implemented usually via continuation methods 
[Kuznetsov, 1995]. In this study, the software 
AUTO [Doedel et al., 2006] has been used for the 
bifurcation analysis along the parameter line p u - 

3. APPLICATION TO A MAPK SIGNALING 
MODULE 

3.1 Model description 

The method presented in Section 2 has been 
applied to an ODE model of a mitogen activated 
protein kinase (MAPK) signaling module. MAPK 
signaling is a recurring motif in cellular signaling 
pathways, and typically appears in a cascade 
involving three levels [Pearson et al., 2001]. 

For this study, we consider a mathematical model 
for the Ras/Raf signaling pathway similar to the 
one presented by Kholodenko [2000]. The inhibi- 
tion of the upstream molecule SOS by activated 
MAPK, the lowest level in the cascade, constitutes 
a negative feedback circuit around the cascade. 
Via the loop-breaking approach, the influence of 
this feedback connection on existence of sustained 
oscillations in kinase activity is analysed. 
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Fig. 1. Illustration of the MAPK cascade model. 
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Table 1. Reaction rates in the MAPK 
cascade model 



The structure of the model is illustrated in Fig. 1. 
The reaction rates as labeled in the figure are 
displayed in Table 1. The concentrations have 
been denoted as xu = [MAPKKK*], x 2 \ = 
[MAPKK*], x 22 = [MAPKK**], x 31 = [MAPK*] 
and X32 = [MAPK**]. The concentrations of un- 
phosphorylated kinases can be computed by con- 
servation laws and the three parameters x\ t , x 2 t 
and x 3t for the total concentrations of the three 
kinases. The difference to the model from Kholo- 
denko [2000] is that the phosphorylation reactions 
3, 4, 7 and 8 are assumed to follow mass ac- 
tion rather than Michaclis-Menten kinetics. This 
is reasonable since the Michaelis-Menten kinetics 
assumes low enzyme concentration compared to 
the substrate, whereas the concentrations of the 
kinases are in a comparable range here. Nominal 
parameter values have been adopted from Kholo- 
denko [2000], and are shown in Table 2 as p n . 

Using the reaction rates from Table 1, the model 
can be written as a system of five ODEs with 20 
parameters: 

in =vi-v 2 

±21 = "3 + V 5 -V4-V 6 

±22 = V4 ~ V 5 (10) 

±31 = V 7 + Vg - V 8 - V 10 

±32 = V s - Vg 
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Fig. 2. Convergence to steady state for parameters 
Po (dashed line) and sustained oscillations for 
parameters pi (solid line). The oscillations 
coexist with an unstable equilibrium (dotted 
line). 

For the nominal parameters po, the model has a 
stable equilibrium x(p n ). Solutions of the model 
converge quickly to the steady state, as depicted 
in Fig. 2. 
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Table 2. Reference parameters p and 
parameters for instability p\ in the 
MAPK cascade model. 



3.2 Parameters for a change in stability properties 

This section describes the application of the 
method presented in Section 2 to the problem of 
finding destabilizing parameters for the MAPK 
cascade model (10). 

The first step is to choose a suitable loop breaking. 
For the MAPK cascade, an intuitive approach is 
to break the loop at the feedback inhibition of 
reaction v\ by MAPK**. Thus we choose h(x) = 
x 32 , to select [MAPK**] as an output, and replace 
X32 by the input u in the reaction rate v\ to obtain 
the dynamics of the open loop system f(x,u,p). 

It can be shown that there is a unique equilib- 
rium of (10) for any parameters in the biologi- 
cally meaningful range. The equilibrium can easily 
be computed numerically. A linearisation of the 
open loop system around this equilibrium point 
and a Laplace transformation gives the transfer 
function G(p, s), whose graph is shown in Fig. 3. 
The set of critical frequencies is minimal with 
a = 3, which can be seen from Fig. 3 by the 
observation that the graph of G(po,juj) encircles 
the origin monotonically. The only positive critical 
frequency is uj c (po) — 0.017s -1 , and we will con- 
sider this frequency in the search for destabilizing 
parameters. The corresponding transfer function 
value is G(po, ju) c (Po)) = 0.12, corresponding to 
the equilibrium x(p ) being stable in the closed 
loop system. 

For the computational approach described in Sec- 
tion 2.5, we chose 7 = 1.5, such that the value of 
the transfer function would have to pass the point 
1 when going from its inital value of 0.12 to 7. The 
optimization method converges to the parameters 
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Fig. 3. Nyquist plots of open-loop MAPK model 
for parameters po (dashed line) and pi (solid 
line). 

pi, which give the desired value G(p\, ju c (p\)) = 
1.5 at a critical frequency oj c (pi) — 0.0065s -1 . 
The parameter values in pi are shown in Table 2. 
The maximal single parameter change from p to 
Pi has been restricted in the numerical implemen- 
tation to be not more than a factor of 5. Even with 
this restriction, parameters leading to sustained 
oscillations have been found. However, 11 out of 
the 20 parameters have been changed by more 
than 20 % to achieve this. 

The graph of G(pi,ju) is shown in Figure 3. For 
the new parameters pi, the graph now encircles 
the point 1. By the argument principle, we see 
that the linearisation of the closed-loop system 
around the equilibrium has some eigenvalues in 
the right half complex plane and is thus unstable. 
The sustained oscillations that appear in this case 
are shown in Fig. 2. 

In conclusion, our method is able to compute 
parameters which render the stable equilibrium 
unstable and thus lead to the emergence of sus- 
tained oscillations. About half of the parameters 
are varied by a non-negligible amount, but all 
variations are within the physiological range. 
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Fig. 4. Bifurcation diagram along the line p^, 
showing stable equilibrium (solid line) , unsta- 
ble equilibrium (dashed line) and amplitude 
of oscillations (circles) . 

3.3 Bifurcation analysis along a line 

Let us now consider the line p^ ~ po + ft{p\ — 
Po). By classical bifurcation analysis with ^ as 
bifurcation parameter, we can see how the system 
changes from the stable to the unstable equilib- 
rium. The resulting bifurcation diagram is shown 
in Figure 4. As expected, there is a Hopf bifurca- 
tion between po and p\, at jj, = 0.664. The evo- 
lution of the limit cycle producing the sustained 
oscillations along the line in parameter space is 
obtained from the bifurcation diagram. 



4. CONCLUSIONS 

We introduced some theoretical tools to inves- 
tigate the existence of parameters for which a 
bifurcation can occur in a dynamical system with 
a feedback circuit. These tools gave rise to a new 
computational method which allows to search for 
parameter values such that the stability proper- 
ties of an equilibrium change in a specific way 
compared to the nominal parameter values. Our 
approach is particularly useful if there are many 
parameters in the system which can be varied si- 
multaneously, and if the contribution of individual 
parameters to stability properties is not obvious. 
The ability to directly handle multiparametric 
variations is a clear advantage compared to using 
only classical bifurcation analysis. 

We have shown the application of the proposed 
method to a model of a MAPK cascade. Using 
relatively small changes to most of the 20 param- 
eters in the model leads to a change from a stable 
equilibrium to an unstable equilibrium with a sta- 
ble limit cycle, producing sustained oscillations. 



Optimal control, stabilization, and nonsmooth analysis, 
pages 135-154. Springer- Verlag, 2004. 

V. Chickarmanc, B. N. Kholodcnko, and H. M. 
Sauro. Oscillatory dynamics arising from compet- 
itive inhibition and multisitc phosphorylation. J. 
Theor. Biol, 244(l):68-76, January 2007. URL 
http : //dx.doi . org/10 . 1016/j . jtbi . 2006 . 05 . 013. 

O. Cinquin and J. Demongcot. Positive and neg- 
ative feedback: striking a balance between neces- 
sary antagonists. J. Theor. Biol, 216(2):229-241, 
May 2002. doi: 10.1006/jtbi.2002.2544. URL 
http : //dx . doi . org/ 10 . 1006/ j tbi . 2002 . 2544. 

B. F. Dibrov, A. M. Zhabotinsky, and B. N. Kholodenko. 
Dynamic stability of steady states and static stabiliza- 
tion in unbranched metabolic pathways. J. Math. Biol., 
15(l):51-63, 1982. 

E. J. Doedel, R. C. Paffcnroth, A. R. Champneys, T. F. 
Fairgrieve, Y. A. Kuznetsov, B. E. Oldeman, B. Sand- 
stede, and X. Wang. AUTO 2000: continuation and 
bifurcation software for ordinary differential equations. 
Concordia University, Montreal, Canada, 2006. 

T. Eifiing, S. Waldherr, F. Allgower, P. Scheurich, 
and E. Bullingcr. Steady state and (bi-) stability 
evaluation of simple protease signalling networks. 
BioSystems, Epub ahead of print, 2007. URL 
http : //dx . doi . org/ 10 . 1016/j . biosystems . 2007 . 01 . 003. 

M. Kaufman and R. Thomas. Emergence of complex 
behaviour from simple circuit structures. Comptes rend, 
biol., 326:205-214, 2003. 

B. N. Kholodcnko. Negative feedback and ultrasensitivity 
can bring about oscillations in the mitogen-activated 
protein kinase cascades. Eur. J. Biochem., 267(6): 1583- 
88, Mar 2000. 

J. Kim, D. G. Bates, I. Postlethwaite, L. Ma, and P. A. 
Iglcsias. Robustness analysis of biochemical network 
models. IEE Proc. Syst. Biol., 153(3):96-104, May 2006. 

Y. A. Kuznetsov. Elements of Applied Bifurcation Theory. 
Springer- Verlag New York, 1995. 

N. I. Markcvich, J. B. Hoek, and B. N. Kholodcnko. Signal- 
ing switches and bistability arising from multisite phos- 
phorylation in protein kinase cascades. J. Cell Biol., 
164(3):353-359, Feb 2004. doi: 10.1083/jcb.200308060. 
URL http : //dx . doi . org/ 10 . 1083/ j cb . 200308060. 

Optimization toolbox. For use with Matlab. The Math- 
Works Inc., 2006. 

G. Pearson, F. Robinson, T. B. Gibson, B. E. Xu, 
M. Karandikar, K. Berman, and M. H. Cobb. Mitogen- 
activated protein (MAP) kinase pathways: regulation 
and physiological functions. Endocr. Rev., 22(2):153- 
183, Apr 2001. 

C. D. Thron. The secant condition for instability in bio- 
chemical feedback-control. 1. The role of coopcrativity 
and saturability. Bull. Math. Biol., 53(3): 383-401, 1991. 

J. J. Tyson and H. G. Othmer. The dynamics of feedback 
control circuits in biochemical pathways. Progr. Theor. 
Biol., 5:2-62, 1978. 

N. Yildirim and M. C. Mackey. Feedback regulation in 
the lactose opcron: a mathematical modeling study and 
comparison with experimental data. Biophys. J., 84(5): 
2841-51, May 2003. 



REFERENCES 

D. Angcli and E. D. Sontag. Interconnections of mono- 
tone systems with steady-state characteristics. In 
M. de Queiroz, M. Malisoff, and P. Wolcnski, editors, 



